% This script starts estimation of the Macro finance model using the SR approach
% By Martin M. Andreasen, August 2019
addpath(genpath(pwd))
close all
clear
%clc
%% User settings
dataOption        = 2;          %1 for using ygap = UGAP,
                                %2 for using ygap = potential output in GDP
                                %3 for using ygap = Cooper and Priestley gap
startYear         = 1961;       %startYear for the estimation. Data available from 1961
modelNum          = 1;          %1 for QTSM without interest rate smoothing
                                %2 for QTSM wiht constant interest rate smoothing

% For inference
theta12OptimalOn  = 0;          %1 for setting LAMBDA in step 3 in an optimal manner. 0 for using the estimates from step 2
AVarTheta1On      = 0;          %1 for computing AVar for theta1 at step 1 and theta11 at step 3, else 0

% The optimizer:
optimizerStep1     = 0;         % optim = 0 for no optimization
optimizerStep3     = 0;         % optim = 1 for the CMAES
                                % optim = 2 for a gradient based optimizer (the LM version)
                                % optim = 3 for the Nelder-Mead.
                                % optim = 4, for first a gradient optimizer and the the Nelder-Mead
optim.numOptim    = 1;          % Number of optimizations - if optimizer = 3 or 4.
optim.sigma       = 0.010;       % step length for CMAES
optim.sigmaMax    = 0.050;       % max length for CMAES
optim.TolFun      = 1e-5;
optim.TolX        = 1e-5;
optim.MaxIter     = 4000;
optim.MaxFunEvals = 4000;


epsValue          = 1D-5;       % Stepsize for the standard errors of Theta1
AVarAdjThetaOn    = 0;          % 1 to adjust factor variance for Theta1, else 0

h0RegimeOn        = 1;          %1 for switching in h0 under P, else 0
hxRegimeOn        = 1;          %1 for switching in hx under P, else 0

%% Default setting
numCPUs           = feature('numcores'); % Number of CPUs
fullPHIon         = 1;          %1 for a full phi, else 0 for diagonal phi
nm                = 2;          %Number of macro factors
nn                = 0;          %Number of latent factors factors
loadDefaultSetting
%factorSignResOn = 2
factorSignResOn = 2

% Loading the data:
loadData = getData(dataOption,startYear);

%% Starting values for the estimated parameters
theta1 = initializeTheta1(fullPHIon,nm,nn);

% Impose restrictions on phi from ReStud paper
theta1 = rmfield(theta1,'phi14');
theta1 = rmfield(theta1,'phi23');
theta1 = rmfield(theta1,'phi32');
theta1 = rmfield(theta1,'phi34');
theta1 = rmfield(theta1,'phi41');
theta1 = rmfield(theta1,'phi43');

if modelNum == 1
    if dataOption == 1
        theta1 = rmfield(theta1,'rhor');
        theta1ValuesStep1 = [   0.000007697957    0.088798954825    0.000001448479    0.006434745357    0.025271178113  ...
            0.003801472813    0.004258045841   -0.485207099637   -0.491153447384   -0.021787206525  ...
            -0.000928735436    0.002855108937   -0.024156823397    0.000103525349    0.000830331327  ...
            0.023919846878   -0.087674048434    0.049990615980    0.018135323127   -0.000302856709 ...
            0.039970173130    0.002512837718   -0.004243300828    0.002117672196    0.000003903177  ]';
        % Q = 0.069862979582060, done by CMAES
        
        % For step 3
        if h0RegimeOn == 0 && hxRegimeOn == 0
            theta1ValuesStep3 = [   -0.000054267465    0.062141749544    0.000001029203    0.020190705171    0.023033327077  ...
                0.002630766803    0.004823063396   -0.252093024059   -0.361085717407    0.015628416992  ...
                -0.000725748632    0.002047542392   -0.015202204828    0.004076399520    0.000672800837  ...
                0.001655995736    0.000039057451    0.001943359316   -0.006279208110   -0.001922424891  ...
                0.011801466352   -0.000719350910   -0.000207358753    0.000926326930    0.000712422918 ]';
            % Q =   0.073725835300225, done by CMAES
        elseif h0RegimeOn == 0 && hxRegimeOn == 1
            theta1ValuesStep3 = [   -0.000052529151    0.061893644107    0.000001015031    0.020127398256    0.023015775710  ...
                0.002649581929    0.004818868026   -0.249851076249   -0.362188831135    0.015830591983  ...
                -0.000727219081    0.002044829664   -0.015189659630    0.004084626447    0.000671770332  ...
                0.001634667604    0.000033630677    0.001939529654   -0.006542854015   -0.001840507312 ...
                0.011472572799   -0.000733829336   -0.000202007946    0.000909291558    0.000711139421 ]';
            % Q= 0.073727371238604, done
        end
    elseif dataOption == 2
        theta1 = rmfield(theta1,'rhor');
        theta1ValuesStep1 = [  0.000727510438    0.561228836590    1.309973087385    0.005818276338    0.033643749878  ...
            0.006207972845   -0.008295617748   -0.048837481175   -9.989960502617   -0.442816157073  ...
            -0.002476332242    0.063225151396   -0.186744386352   -0.051991162114   -0.001266669316  ...
            0.055013577028    0.190841576621    0.160741443933   -0.030392501464   -0.021191895909  ...
            0.007433157416    0.001924109609   -0.006640332880    0.003645634498    0.000124700725    ]';
        % Q = 0.076035080792074, done
        
        % For step 3
        if h0RegimeOn == 0 && hxRegimeOn == 0
            theta1ValuesStep3 = [    -0.000334635113    0.148027313324    1.744012928324    0.009536726439    0.033435319775  ...
                0.032824466428   -0.007925034068    0.355601071136   -9.928400276365   -0.095084756859  ...
                -0.000054338968    0.016683195355    0.245461180748   -0.006481673078    0.000399557182  ...
                0.001608347989   -0.002013128380    0.062737252892   -0.001162363419    0.000958539536 ...
                0.003203224780   -0.000029158204    0.000122453277    0.001232349536    0.000651932900 ]';
            % Q =   0.083714083339630, done
        elseif h0RegimeOn == 0 && hxRegimeOn == 1
            theta1ValuesStep3 = [     -0.000898461204    0.250851649778    1.738812153685    0.000082749752    0.026741235717  ...
                0.027608889948    0.002757913198    0.460302792874   -9.897605283470   -0.173810842545  ...
                -0.000197058907    0.025329589583    0.328941504271   -0.013121404295    0.000231383379  ...
                0.001593344261   -0.005290459065    0.054468867092   -0.001185198844    0.000771546587  ...
                0.003115767926   -0.000029576045    0.000076376834    0.001193961938    0.000650259382  ]';
            % QStep3 =  0.084088556155478
        end        
    end
end
theta1           = values2struct(theta1ValuesStep1,fieldnames(theta1));
theta1StartStep3 = values2struct(theta1ValuesStep3,fieldnames(theta1));

% Getting bounds and Insigma
[upperBounds,lowerBounds,Insigma] = theta1BoundsInsigma(nm,nn);

% Defining the struct with the setup specification for step 1
constructSetupStep1

%% Computing the objective function theta1
tic
[Q,output,errorMes] = objectFunc(struc2values(theta1),setupStep1);
toc
%% Step 1 of the SR approach
% The optimization
optim.optimizer = optimizerStep1;
[QStep1,theta1Step1,outputStep1] = step1SR(theta1,setupStep1,optim);
theta1 = theta1Step1;

% The standard errors
AVarTheta1Step1  = getAVarTheta1andFactors(theta1Step1,setupStep1,epsValue);
%plotFactors(setupStep1,outputStep1,AVarTheta1Step1)
%% Step 2 of the SR approach
[theta2Step2,AVarTheta2Step2] = step2SR_threshold_macro(AVarTheta1Step1,outputStep1,bandwidthTstep2,setupStep1.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
%mprStep2 = getMPR(outputStep1,theta2Step2,AVarTheta2Step2,nm*2+nn,h0RegimeOn,hxRegimeOn);
%% Step 3 of the SR approach
namesTheta12 = cell((2*nm+nn)*(2*nm+nn+1)/2,1);
idx = 0;
for i=1:2*nm+nn
    for j=1:i
        idx = idx + 1;
        namesTheta12{idx} = ['sigma',num2str(i),num2str(j)];
    end
end
[theta12Step3,AVarTheta12Step3] = step3SR_theta12(namesTheta12,theta1Step1,AVarTheta1Step1,theta2Step2,AVarTheta2Step2,theta12OptimalOn);

% Estimation of theta11 - using Robust version
if exist('theta1StartStep3','var')
    theta1Start = theta1StartStep3;
    if AVarTheta1On  == 1
        disp('for standard errors')
        theta1Start = rmfield(theta1Start,'phi22');
        setupStep1.calibrateTheta1.phi22 = theta1StartStep3.phi22;
        setupStep1.selectTheta1          = fieldnames(theta1Start);
        setupStep1.upperBoundsValues    = struc2values(upperBounds,setupStep1.selectTheta1);
        setupStep1.lowerBoundsValues    = struc2values(lowerBounds,setupStep1.selectTheta1);
        setupStep1.InsigmaValues        = struc2values(Insigma,setupStep1.selectTheta1);
    end
else
    theta1Start = theta1Step1;
end
optim.optimizer = optimizerStep3;
[QStep3,theta11Step3,outputStep3,setupStep3] = step3SR_theta11(theta12Step3.Robust,namesTheta12,theta1Start,...
    setupStep1,upperBounds,lowerBounds,Insigma,optim);

% Standard errors of theta11: The asymptotic standard errors
AVarTheta11Step3 = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step2,AVarTheta12Step3,...
    setupStep1,setupStep3,epsValue/10,theta12OptimalOn);

% Re-estimating P dynamics
[theta2Step3,AVarTheta2Step3] =step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);

% Enforce consistency of sigma between Q and P
posTheta12Step2 = zeros(1,length(namesTheta12));
for i=1:length(namesTheta12)
    theta2Step3.(namesTheta12{i}) = theta2Step2.(namesTheta12{i});
    posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step3.names,namesTheta12(i)));
end
AVarTheta2Step3.VarRobust(posTheta12Step2,posTheta12Step2) = AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2);

mprStep3 = getMPR(outputStep3,theta2Step3,AVarTheta2Step3,nm*2+nn,h0RegimeOn,hxRegimeOn);
testRegimeStep3 = waldTestRegimeShift(theta2Step3,AVarTheta2Step3,2*nm+nn);

%% Do the yield curve decomposition
if h0RegimeOn == 0 && hxRegimeOn == 0
    yieldTP = yieldCurveDecom(setupStep3,outputStep3,theta2Step3);
else
    numSimTP = 1D4;
    yieldTP  = yieldCurveDecom_threshold(setupStep3,outputStep3,theta2Step3,numSimTP);
end

% mean(yieldTP.termPremia(end,:))
% std(yieldTP.termPremia(end,:))
%% Simulating the model
thresholdLevel          = AVarTheta2Step3.thresholdLevel;
momentsData             = getMomentsData(loadData,outputStep3,thresholdLevel);
momentsModel            = simulator_threshold(theta2Step3,outputStep3,thresholdLevel,numSim);
momentsModel_finite     = simulator_threshold_finite(theta2Step3,outputStep3,thresholdLevel,1000);
expectedExcessReturn.h1 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,1);
expectedExcessReturn.h3 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,3);

%% Do plots
plotFactors(setupStep3,outputStep3,AVarTheta11Step3)
plotFit(setupStep3,outputStep3)
plotUnconditionalMoments(momentsModel,momentsData)
plotForecastRegressions(momentsModel,momentsData)
plotYieldCurveDecom(yieldTP,setupStep3)

%% Displaying the results
outTheta1          = displayResults(theta1Step1,AVarTheta1Step1.VarRobust);
outTheta11         = displayResults(theta11Step3,AVarTheta11Step3.VarRobust);
theta2Step3Reduced = reduceTheta2(theta2Step3,AVarTheta2Step3.names);
outTheta2Reduced   = displayResults(theta2Step3Reduced,AVarTheta2Step3.VarRobust);

disp('optimum at step 1')
printToScreen(struct2array(theta1Step1));
disp('optimum at step 3')
theta1Step3Start  = theta11Step3;
for i = 1:size(namesTheta12,1)
    theta1Step3Start.(namesTheta12{i}) = setupStep3.calibrateTheta1.(namesTheta12{i});
end
printToScreen(struct2array(theta1Step3Start));

%% Saving the results in a mat file
if modelNum == 1
    modelName = 'QTSM_nn0_ReStudPhiQ_outputGap_';
elseif modelNum == 2
    modelName = 'QTSMr_nn0_ReStudPhiQ_outputGap_';
end
filename = [modelName,'dataOption',num2str(dataOption),num2str(startYear),'_SE',num2str(epsValue),'_h0RegimeOn',num2str(h0RegimeOn),'_','hxRegimeOn',num2str(hxRegimeOn),...
    '_numBootStep2_',num2str(numBootStep2)];
save([pwd,'\Results\',filename,'.mat'],'QStep1','theta1Step1','outputStep1','AVarTheta1Step1',...
    'theta2Step2','AVarTheta2Step2',...
    'theta11Step3','outputStep3','setupStep3','AVarTheta11Step3','theta2Step3','AVarTheta2Step3',...
    'yieldTP','momentsModel','momentsModel_finite','momentsData','mprStep3','testRegimeStep3','expectedExcessReturn')


